Primordial clouds under the influence of background radiation field 
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Our goal is to study the effects of the UV radiation from the first 
stars, quasars and hypothetical Super Heavy Dark Matter (SHDM) parti- 
cle decays on the formation of primordial bound objects in the Universe. 
We trace the evolution of a spherically symmetric density perturbation in 
the Lambda Cold Dark Matter and MOND model, solving the frequency- 
dependent radiative transfer equation, non-equilibrium chemistry, and one- 
dimensional gas hydrodynamics. We concentrate on the destruction and 
formation processes of the Hi molecule, which is the main coolant in the 
primordial objects. 

PACS numbers: 95.30.Lz, 95.30.Sf, 98.35.Mp, 98.80.Bp 



1. Introduction 

Formation of primordial objects such as young galaxies and globular 
clusters, in the early universe is a fundamental problem in the modern cos- 
mology. There was a rapid progress of observations in the last years. It 
gave theoreticians invaluable possibility to compare their theories with the 
observations. Now, we can say that the framework of the structure forma- 
tion is well known. However, there are still many unanswered questions, 
especially regarding the influence of the background radiation field from the 
first objects to the star formation rate. In our work we are trying to answer 
this question. 
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The existence of the first objects is a direct consequence of the growth 
of the primordial density fluctuations. At the beginning, there are linear 
density perturbations which expand with the overall Hubble flow. Subse- 
quently, these perturbations can grow and form primordial clouds. Clouds 
with enough density contrast decouple from this flow and start to collapse. 
The kinetic energy of the infalling gas is dissipated through shocks and the 
cloud becomes pressure supported. The further evolution of the cloud is 
determined by its ability to cool sufficiently fast. Clouds which could not 
cool fast enough will stay in a pressure-supported stage and will not form 
any stars. The existence of the efficient cooling mechanism is necessary to 
continue the collapse of the cloud, its subsequent fragmentation and star 
formation. 

In our work we are interested in the first generation of stars, say III 
population, which are still forming in the low mass clouds when first lumi- 
nous objects already exist. These objects are made from the primordial gas 
so they are metals free. It is simply because the first stars did not have 
much time to produce them. These objects could be irradiated by the UV 
and X-rays radiation produced by the first stars, quasars and hypothetical 
Super Heavy Dark Matter (SHDM) particle decays (p], [2]). 

In the absence of metals, the most important cooling mechanism for low- 
mass primordial clouds is so called 'i/2 cooling', i.e. cooling by radiation 
of excited rotational and vibrational states of H2 molecule. We can look at 
this mechanism as the collisional excitation of the H2 molecule, subsequent 
spontaneous de-excitation and photon emission. Emitted photon can further 
escape from the cloud and take away kinetic energy of the collapsing cloud. 
The presence of initial mass fraction of the molecular hydrogen H2 of only 
10~ 6 is enough to trigger the final collapse of low mass clouds. 

Molecular hydrogen is fragile and can easily be photodissociated by pho- 
tons with energies of 11.26 - 13.6 eV (Lyman and Werner bands)([3j) through 
the following reactions 

H 2 + his -> HI -» 2H. 

Destruction of the H2 would stop collapsing of the low mass clouds and 
decrease star formation rate. 

In primordial gas cloud, H2 molecules can form mainly through the 
coupled reactions 

e" + H -> H~ + hv 
H- + H — H 2 + e~, 

in which the electrons act only as a catalyst. X-rays radiation can increase 
production of the H2 by enhancement of free electrons fraction. 
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From above we see that the UV and X-rays radiation background alters 
the subsequent growth of cosmic structures. It regulates star formation rate 
so it has important implication to the re-ionization history of the Universe 
([!]). In addition photoionization by UV photons can be regarded as the 
process which lead to the inhibition of the low mass galaxies formation, so 
it could explain the so called 'low-mass galaxies overproduction problem' 
in the hierarchical bottom-up theory of galaxy formation. It is therefore 
crucial to determine quantitatively the consequences of radiation feedback 
on the formation of early generation objects. 

Feedback of the UV background to the collapse of the spherically sym- 
metric primordial gas cloud in CDM model was studied by many authors 
(0, [S], [7], [8], [9]). However, their calculations were simplified. They have 
been using self-shielding function ([10]) instead of solving the full transfer 
equation. Our approach includes solving the frequency-dependent radiative 
transfer equation in the 'exact' way. We try to check the evolution of the 
collapsing cloud in the ACDM and MOND model. 



2. Basic equations 

In Cold Dark Matter models and in case of spherical symmetry, the 
collapse may be described by the following equations: 



dM 

dr 
dr 

~dl 



Airr 2 g, (1) 

(2) 



2,.+- W 



dv 9 dp GM(r 

— = -^r 2 — 

dt dM r 

du p dg A 
dt q 2 dt g 

where r is the radius of a sphere of mass M, u is the internal energy per 
unit mass, p is pressure and g is mass density. Here, Eq.([T]) is the continuity 
equation, ([2]) and (|3|) give acceleration and (jU) accounts for energy conser- 
vation. The last term in Eq.dH) describes gas cooling/heating, where A is 
energy absorption (emission) rate per unit volume. 

A consist of two parts, that is, the chemical cooling A c hem and the ra- 
diative cooling A ra( i: 

A (r) = A chem (r) + A rad (r) . (5) 
The A c /j em can be written as 

Achem = Q (O) 
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where e c hem is the chemical binding energy per unit mass. The radiative 
cooling A ra( i (r) is obtained by solving radiative transfer equation. 
We use the equation of state of perfect gas 

p =(j-i)gu, (7) 

where 7 = 5/3, as the primordial baryonic matter after recombination is 
assumed to be composed mainly of monoatomic hydrogen and helium, with 
the fraction of molecular hydrogen H2 always less than 10~ 3 . 
In the case of modified gravity [TT], equation © will change: 

dv . 2 dp ( GMjr) g H \ fa . 

T77 = ~ 47rr 777,7 ~ 9H ~ aof [ s , (8) 

at dM \ aor z ao J 

where f(x) is some function, asymptotically equal to x for x 3> ao and to 
yja$x for x <C ao, while gn may be expressed as 

9h = \h 2 [(z + ifn b + 2 ((* + l) 4 O r - n A )] r. (9) 

Here, Hq is the current value of the Hubble parameter, fij,, f2 r and £l\ are 
the current fractions of baryons, radiation and dark energy in terms of the 
critical density of the Universe, and z is the redshift. We have made similar 
assumption like in [12] and we apply MOND to net acceleration over Hubble 
flow only. 

3. Radiative transfer equation in spherical symmetry 

The time- independent, non-relativistic equation for radiation transport 
in spherical geometry can be written: 

+ l — (L^ = g { T}v (r i i i )- Xt/ (r)i v y i (10) 
or r o\i 

where I v = l v (r, fi) is the intensity of radiation of frequency u, at radius r 
and in the direction /i = cos#, where 9 is the angle between the outward 
normal and photon direction. r\ v (r, \i) and \v (r) are the total emissivity 
and opacity at frequency v. 

r] v (r, n) and \v ( r ) can De expressed as: 

Vv (r,p) = vl (r) + Vl (r,p) (H) 
Xu{r) = K u {r)+a v {r), (12) 

where rf v , rjf,, k u and o v are thermal emissivity, scattering emissivity, true 
absorption coefficient and scattering coefficient at frequency v respectively. 
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Following Hummer and Rybicki [13J, we introduce a more convenient 
system of coordinates (r, p) rather than (r, fj,) defined by the transformation 
formulae 



(r, fi) — ► \r,p = ryl — fi 2 J for — 1 < < 1. (13) 

For a given radius r the 'impact' parameter p can vary between and r. 
Because the parameter p cannot distinguish between \i > and n < 0, the 
radiation intensity 7j, have to be separated into outward 7+ and inward 7~ 
directed intensity: 

7+ = J„(r,p) /x>0 (14) 
7- = 7,(r,p) /x<0. (15) 



Now, in the new system of coordinates Eq. (|10p can be written as two 
separate equations 

QJ+ 

I l ~q~ = ~eXv {It - S v (r,p)} (16) 
57" 

V-~q~ = eXu{I u -S„(r,p)}, (17) 



where fi = yl — ^ > and S v (r, p) = r\ v (r, p) / \ v (r) is the source func- 
tion. 

If we introduce 1141 



Ju(r,p) = ^{I+(r,p) + I-(r,p)) (18) 
h u (r,p) = ^(I+(r,p)-I-(r,p)), (19) 
than Eqs. (|16p and (|17p will get the form 

» djv £ P) = -QX,K(r,p) (20) 

V — V -q^ — = -QXu{ju(r,p) - S v (r,p)}. (21) 

To solve these equations we need also a boundary condition. We consider 
a spherical envelope with an inner boundary at radius r = r m i n and an outer 
boundary t — ^max- The cloud is immersed in an external time dependent, 
isotropic radiation field. We obtain the following boundary condition 



' milt 



K (r min ,p) = a„(p,t) 0<p< r r , 

K (p,p) = r min <p< r max (22) 

jv { r maxiP) ~ h u {r max , 

p) = Pu(p,t) < p < r r , 



' max ) 



G 



Primordial printed on February 5, 2008 



where a u and (5 V describe the changes of the radiation at frequency u im- 
pinging upon the inner and the outer boundary with the time t. 

By the standard procedure one obtains from Eq. (|10p the zeroth and the 
first moment equations [15] 

d(f u (r)JAr)) , 3/„ (r) - 1 

0^ 1 J v [r) + QX v H v {r) = (23) 

dH u (r) 2H u (r) t 

x ' + — ^ + QKuJu (r) - evt = 0, 24 
or r 

where f v (r) = K v (r) / J v (r) is the Eddington factor. J u , H v , K v are the 
zeroth, first and second moment of the radiation field at frequency v: 

Ju{r) = [ j„dfi (25) 



o 

H v {r) = { Kfidfi (26) 
Jo 

K u (r) = f j u ^ 2 dfi. (27) 
Jo 

We can cast Eq. (l23p into a more convenient form, by introducing a 
sphericality factor q v , defined in the following way: 



q v (r) = exp 



r / 1 \ dr' 
3 



U{r')J r> 



(28) 



where r c is the core radius, that is, the inner boundary of the medium. 
Using this factor and making some simplification we can rewrite Eqs. (l23p 
and dMI) as 

d{fv{r)qv{r)J u (r)) 

= -pXuqv{r)H u {r) (29) 



dr 

d (H u (r) r 2 ) 



r 2 p (jil - k u J u (r)) . (30) 



dr 

The corresponding boundary conditions for the moment equations (|29p 
and (|30p ) can be written as follows [16] 

H u (r min ) = [ a v pdp (31) 



o 

Ju(r ma x) = / Pvd[i+ / K(r max ,p)dfi. (32) 
Jo Jo 

We see that these boundary condition can be determined only after 
Eqs. lpOj) and (pi]) have been solved for j v and h u . 
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In the following sections we shall refer to Eqs. (l20p and (|21|) along with 
the boundary conditions (|22p as system I equations. The moment equations 
(|29p and (|30p along with the boundary condition (|3ip . (|32|) will henceforth 
be referred as system II equations. 

Solution of the radiative energy transfer equations gives us information 
about the luminosity function L (r): 

L (r) = 16vr Vtf (r) , (33) 

where 

H(r)= / H v (r) di/. (34) 



With the given luminosity we can calculate radiative part of cooling 
function A ra rf: 

1 9L(r) 

*-< r >— s?-sr- (35) 

It can be used further in the hydro dynamical equations. 



4. Chemical reactions and thermal effects 

In our calculations, we include all of the relevant thermal and chemical 
processes in the primordial gas. We have taken into account nine species: 
H, H~, H + , He, He + , He ++ , H2, H^ and e~. Their abundance varies with 
time due to chemical reactions, ionization and dissociation photoprocesses. 
The chemical reactions include such processes as ionization of hydrogen 
and helium by electrons, recombination of ions with electrons, formation of 
negative hydrogen ions, formation of H2 molecules, etc. A full list of the 
relevant chemical reactions and appropriate formulae is given in [17j . 

The time evolution of the number density of component is described 
by the kinetic equation: 

j 9 9 9 

-77 = Y a lmih m nin m + Y bjiUjUj. (36) 
dt 1=1 m=l 3=1 

The first component on the right-hand side of this equation describes the 
chemical reactions, and the other one accounts for photoionization and pho- 
todissociation processes. Coefficients k[ m are reaction rates, quantities K n 
are photoionization or photodissociation rates, and a; m , and bji are numbers 
equal to 0, ±1 or ±2 depending on the reaction. All reaction rates, as well 
as photoionization and photodissociation rates are given in |17j . 

The cooling (or heating) function A includes effects of collisional ion- 
ization of H, He and He + , recombination to H, He and He + , collisional 
excitation of H and He + , Bremsstrahlung, Compton cooling and cooling by 
de-excitation of H2 molecules. The heating/cooling rates are given in |17j . 
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5. Numerical strategy and initial conditions 



In the simulations we used the code described in [17], based on those 
presented by [18] and [19]. This is a standard, one-dimensional, second- 
order accurate, Lagrangian finite-difference scheme, slightly modified to our 
purposes. Our code supports flat and non-flat CDM and ACDM models as 
well as Milgrom's Modified Newtonian Dynamics. In the second case it was 
necessary to do significant changes in the initial conditions. We starte our 
calculations at the end of the radiation-dominated era instead of z = 500. 
For Qf, = Qm = 0.02/ h 2 , z e q — 485 as given by the formula provided by [20J , 
z eq = 2.50 x W A Q h 2 @^, where 9 2 . 7 = T 7 /2.7 K, assuming h = 0.72 and 
T 7 =2.7277 K. We assume that in MOND, like in the standard cosmology, 
initial overdensities may grow only in the matter-dominated era. In both 
cases we use our own code to calculate the initial chemical composition and 
initial gas temperature. Initial overdensities may be calculated from the 
matter power spectrum which may be obtained e.g. using the CMBFAST 
program by [21j. 

We apply the initial density profiles in the form of a single spherical 
Fourier mode, also used by [19] 



where g c is the critical density of the Universe, g c = 3H 2 /8ttG, with H 
being the actual value of the Hubble parameter. 

For this profile, we can distinguish two radius values, Rq and R z , which 
correspond to the first zero and the first minimum of the function sin(/cr) /kr, 
respectively. Inside the sphere of radius i?o = ir/k which contains mass Mq, 
the local density contrast is positive. The mass Mq and the radius Rq will be 
referred to as the cloud mass and the cloud radius, respectively. The local 
density contrast is negative for R z > r > Rq, with average density contrast 
vanishing for the sphere of radius R z = 4. 49341/ A; and mass M z . The shell 
of radius R z will expand together with the Hubble flow, not undergoing 
any additional deceleration. This is why this profile is very convenient in 
numerical simulations. It eliminates the numerical edge effects and the 
mentioned shell simply follows the Hubble expansion of the universe. Thus, 
it can be regarded as the perturbation boundary, with mass M z referred to 
as the bound mass. 

It is worth to note that for radii not greater than 3/4Rq this profile is 
very similar to the Gaussian profile 



g b (r) = O fe £ c (l + 5 



sin kr 



(37) 



At 




(38) 



Primordial printed on February 5, 2008 



9 



with Rf = 1/2R . 

As the initial velocity, we use the Hubble velocity: 

v{r) = Hr. (39) 

Our numerical computational procedure to trace the dynamical evolu- 
tion of the primordial cloud under the UV background is following: 

• We solve the hydrodynamic equations of motion along with equations 
for energy conservation, ionization, and dissociation of molecular and 
atomic species. From the solution of these equations we can calculate 
the properties and distribution of the absorbing components k u (r), 
a v (r) and the thermal emissivity rjl(r). It will give us the initial 
value of the source function S u . 

• With the initial value of S v , we solve the system I equations for the 
geometrical factors f v and q u . 

• We solve the system II equations for J u and H v . 

• The source function is independent from the radiation field (J„) only 
when there is no light scattering. So, in general we have to update S u 
and solve the system I once more. 

• Iterative procedure between system I and II is continued until conver- 
gence. 

• We calculate luminosity L (r) from Eq. (|33p and than cooling function. 

• We update abundance of different species and number densities of each 
atomic and molecular state. 

• We repeat all of the above steps. 

6. Results 

Detailed results of our simulations will be described in a paper, which 
is in preparation. 
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